/*
 * Translation to C from nbody0 example in Sverre Aarseth's website.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define min(x,y) (x<y ? x : y)

float *X, *XDOT;

float eta,dtnext,tcrit,eps2;
long Np;

void allocate_arrays( long N )
{
        X     = (float *)malloc( sizeof(float)*4 * N );
        XDOT  = (float *)malloc( sizeof(float)*4 * N );
}

void read_ic( )
{
        FILE *f;
        f = fopen( "../IC/init.txt", "r" );
    
        fscanf( f, "%ld%f%f%f%f\n", &Np, &eta, &dtnext, &tcrit, &eps2 );
        
        allocate_arrays( Np );
        
        for( int j=0; j<Np; j++ )
        {       
                fscanf( f, "%f", &X[4*j+3] );
                
                for( int k=0; k<3; k++ )
                        fscanf( f, "%f", &X[4*j+k] );
        
                for( int k=0; k<2; k++ )
                        fscanf( f, "%f", &XDOT[4*j+k] );
        
                fscanf( f, "%f\n", &XDOT[4*j+2]);
        
        }
        
        fclose( f );

}

void dump_phasecoords( FILE *phstream )
{
        for( long j=0; j<Np; j++ )
        {
                fprintf( phstream, "%f ", X[4*j+3] );
                
                for( int k=0; k<3; k++ )
                        fprintf( phstream, "%f ",  X[4*j+k] );
                
                for( int k=0; k<2; k++ )
                        fprintf( phstream, "%f ", XDOT[4*j+k] );
                
                fprintf( phstream, "%f\n", XDOT[4*j+2]);
        }
}

void energy( float *eK, float *uK )
{
        float K = 0.0,U = 0.0;
        
        for( long i=0; i<Np; i++ )
        {
                K += X[4*i+3]*(pow(XDOT[4*i+0],2)+pow(XDOT[4*i+1],2)+pow(XDOT[4*i+2],2));
                
                for( long j=i+1; j<Np; j++ )
                        U += X[4*i+3]*X[4*j+3]/sqrt( pow(X[4*j+0]-X[4*i+0],2) +
                                                     pow(X[4*j+1]-X[4*i+1],2) + 
                                                     pow(X[4*j+2]-X[4*i+2],2) + eps2 ); 
        }
        
        *eK = 0.5*K;
        *uK =     U;
}

void scale_units( )
{
    float T, U;
    
    float M = 0.0;
    
    long i;
    for( i=0; i<Np; i++ )
        M += X[4*i+3];
    
    for( i=0; i<Np; i++ )
        X[4*i+3] /= M;
    
    
    energy( &T, &U );
    
    float q = sqrt( 0.5*U/T );
    
    float E0 = -1.0/4.0;
    
    float S  = E0/( q*q*T - U );
    
    int k;
    for( i=0; i<Np; i++ )
    for( k=0; k< 3; k++ )
    {
        X[4*i+k] *= 1.0/sqrt(S);
        
        XDOT[4*i+k] *= q*sqrt(S);
    }
    
    M = 0.0;
    for( i=0; i<Np; i++ )
        M += X[4*i+3];

    energy( &T, &U );
    
    fprintf( stderr, "%2.4f  %2.4f\n", q,   S   );
    fprintf( stderr, "%2.4f  %2.4f\n", M, T - U );
    
    printf( "%ld%f%f%f%f\n", Np, eta, dtnext, tcrit, eps2 );
    
    dump_phasecoords();
    
}
    
int main( void )
{
    read_ic();
    
    scale_units();
    
    return 0;
}
